{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Monte Carlo Method\n",
    "\n",
    "As a simple example of a Monte Carlo method, we will approximate the value of $\\pi$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pyopencl as cl\n",
    "import pyopencl.array\n",
    "import pyopencl.clrandom"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "ctx = cl.create_some_context()\n",
    "queue = cl.CommandQueue(ctx)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Boilerplate for Random Number Generator"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "generator_preamble = \"\"\"\n",
    "#include <pyopencl-random123/philox.cl>\n",
    "\n",
    "typedef union {\n",
    "    uint4 v;\n",
    "    philox4x32_ctr_t c;\n",
    "} philox4x32_ctr_vec_union;\n",
    "\n",
    "\n",
    "uint4 philox4x32_bump(uint4 ctr)\n",
    "{\n",
    "    if (++ctr.x == 0)\n",
    "        if (++ctr.y == 0)\n",
    "            ++ctr.z;\n",
    "    return ctr;\n",
    "}\n",
    "\n",
    "uint4 philox4x32_gen(\n",
    "        uint4 ctr,\n",
    "        uint2 key,\n",
    "        uint4 *new_ctr)\n",
    "{\n",
    "    philox4x32_ctr_vec_union result;\n",
    "    result.c = philox4x32(\n",
    "        *(philox4x32_ctr_t *) &ctr,\n",
    "        *(philox4x32_key_t *) &key);\n",
    "    *new_ctr = philox4x32_bump(ctr);\n",
    "    return result.v;\n",
    "}\n",
    "\n",
    "float4 philox4x32_f32(\n",
    "        uint4 ctr,\n",
    "        uint2 key,\n",
    "        uint4 *new_ctr)\n",
    "{\n",
    "    *new_ctr = ctr;\n",
    "    return\n",
    "        convert_float4(philox4x32_gen(*new_ctr, key, new_ctr))\n",
    "        * 2.3283064365386963e-10f;\n",
    "}\n",
    "\"\"\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Reduction Code"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Complete the sampler code:\n",
    "\n",
    "```\n",
    "mc_preamble_src = \"\"\"\n",
    "\n",
    "#include <pyopencl-complex.h>\n",
    "\n",
    "float compute_sample(int i, unsigned int k1)\n",
    "{\n",
    "    uint4 ctr = { 0, 1, 2, 3 };\n",
    "    uint2 key2 = { i, k1 };\n",
    "    float4 rng_res = philox4x32_f32(ctr, key2, &(ctr));\n",
    "    ...\n",
    "}\n",
    "\"\"\"\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#clear\n",
    "mc_preamble_src = \"\"\"\n",
    "\n",
    "#include <pyopencl-complex.h>\n",
    "\n",
    "float compute_sample(int i, unsigned int k1)\n",
    "{\n",
    "    uint4 ctr = { 0, 1, 2, 3 };\n",
    "    uint2 key2 = { i, k1 };\n",
    "    float4 rng_res = philox4x32_f32(ctr, key2, &(ctr));\n",
    "    \n",
    "    cfloat_t samp0 = cfloat_new(rng_res.s0, rng_res.s1);\n",
    "    cfloat_t samp1 = cfloat_new(rng_res.s2, rng_res.s3);\n",
    "    \n",
    "    float result = 0;\n",
    "    if (cfloat_abs(samp0) <= 1)\n",
    "        result += 1;\n",
    "    if (cfloat_abs(samp1) <= 1)\n",
    "        result += 1;\n",
    "        \n",
    "    return result;\n",
    "}\n",
    "\"\"\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from pyopencl.reduction import ReductionKernel"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Syntax:\n",
    "\n",
    "`ReductionKernel(context, dtype, netural, reduce_expr, map_expr, arguments)`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#clear\n",
    "rknl = ReductionKernel(ctx, np.float32,\n",
    "        neutral=\"0\",\n",
    "        reduce_expr=\"a+b\", map_expr=\"compute_sample(i, k1)\",\n",
    "        arguments=\"unsigned int k1\", preamble=generator_preamble+mc_preamble_src)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "3.14154656\n"
     ]
    }
   ],
   "source": [
    "n = 100000000\n",
    "\n",
    "nsamples_accepted = rknl(15, range=slice(n), queue=queue).get()\n",
    "nsamples = 2*n\n",
    "approx_pi = 4 * nsamples_accepted/nsamples\n",
    "\n",
    "print(approx_pi)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.1+"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}
